function acqResults = acquisition_bladerf(longSignal, settings)
%% mod by xiaoyeyimier
% date : 2021-10-12
% vers : 1.0.1
% func : adjust to deal an IQ data for bladerf
% contact : xiaoyeyimier
% ======================================================================= %
%--------------------------------------------------------------------------
%                           SoftGNSS v3.0
% 
% Copyright (C) Darius Plausinaitis and Dennis M. Akos
% Written by Darius Plausinaitis and Dennis M. Akos
% Based on Peter Rinder and Nicolaj Bertelsen
%--------------------------------------------------------------------------
%CVS record:
%$Id: acquisition.m,v 1.1.2.12 2006/08/14 12:08:03 dpl Exp $

%Function performs cold start acquisition on the collected "data". It
%searches for GPS signals of all satellites, which are listed in field
%"acqSatelliteList" in the settings structure. Function saves code phase
%and frequency of the detected signals in the "acqResults" structure.
%
%acqResults = acquisition(longSignal, settings)
%
%   Inputs:
%       longSignal    - 11 ms of raw signal from the front-end 
%       settings      - Receiver settings. Provides information about
%                       sampling and intermediate frequencies and other
%                       parameters including the list of the satellites to
%                       be acquired.
%   Outputs:
%       acqResults    - Function saves code phases and frequencies of the 
%                       detected signals in the "acqResults" structure. The
%                       field "carrFreq" is set to 0 if the signal is not
%                       detected for the given PRN number. 
 

%% Initialization 
% Initialize intermediate variables
fs = settings.samplingFreq;                                                % get sampling frequency
ts = 1 / fs;                                                               % get sampling time 
CodeRate = settings.codeFreqBasis;
CodeLen = settings.codeLength;
BitsPerCode = round(fs /(CodeRate / CodeLen));                             % Find number of samples per spreading code
signal0DC = longSignal - mean(longSignal);                                 % DC filter
phasePoints = (0 : (BitsPerCode-1)) * 2 * pi * ts;                         % Find phase points of the local carrier wave 
numberOfFrqBins = round(settings.acqSearchBand * 2) + 1;                   % Number of the frequency bins for the given acquisition band (500Hz steps)
caCodesTable = makeCaTable(settings);                                      % Generate all C/A codes and sample them according to the sampling freq.
frqBins     = zeros(1, numberOfFrqBins);                                   % Carrier frequencies of the frequency bins
% Initialize acqResults
acqResults.carrFreq     = zeros(1, 32);                                    % Carrier frequencies of detected signals
acqResults.codePhase    = zeros(1, 32);                                    % C/A code phases of detected signals
acqResults.peakMetric   = zeros(1, 32);                                    % Correlation peak ratios of the detected signals
acqResults.acq_Sat_No   = 0;
% disp some info
fprintf('(');                                                              % Two-dimensional correlation peaks were recorded 
fig = figure();                                                            % for debugging and observation
%% Correlate signals   
for PRN = settings.acqSatelliteList                                        % Searches for the specified list of satellites
    axtemp = subplot(4,8,PRN);                                             % Gen a new axes
    hold(axtemp,'On');                                                     % hold on
    title(axtemp,['PRN:',num2str(PRN)]);                                   % name the axes
    view([45,45]);                                                         % Adjust the Two-dimensional correlation peaks to a suitable view
    caCodeFreqDom = conj(fft(caCodesTable(PRN, :)));                       % Perform DFT of C/A code
    results=zeros(numberOfFrqBins,BitsPerCode);                            % Make the correlation for whole frequency band (for all freq. bins)
    for TempF = 1:numberOfFrqBins                                          % Serial search on Doppler Frequency shift
        frqBins(TempF) = settings.IF - ...                                 % Generate carrier wave frequency grid (0.5kHz step)
                               (settings.acqSearchBand/2) * 1000 + ...
                               0.5e3 * (TempF - 1);
        SinC = sin(frqBins(TempF) * phasePoints);                          % Gen Local replica
        CosC = cos(frqBins(TempF) * phasePoints);
        LO = SinC + 1j*CosC;
        for ms_cnt=1:settings.acq_Acc_ms_No                                % Stripping carrier (Incoherent mode)
            signal=longSignal(BitsPerCode*(ms_cnt-1)+1:...                 % read signal per 1 ms
                ms_cnt*BitsPerCode);
            if(settings.OnlyI)                                             % remove carrier (single)
                I1      = SinC .* signal;                                  % process the I branches separately 
                Q1      = CosC .* signal;                                  % process the Q branches separately
                IQfreqDom1 = fft(I1 + 1j*Q1);                              % remove carrier (IQ)
            else
                IQ1      = LO .* signal;
                IQfreqDom1 = fft(IQ1);
            end
            convCodeIQ1 = IQfreqDom1 .* caCodeFreqDom;                     % Multip. in the F domain = Corr in T domain
            acqRes1 = abs(ifft(convCodeIQ1)) .^ 2;                         % IDFT to get correlation results
            results(TempF, :) = results(TempF, :) + acqRes1;
        end
    end 
    [~,frequencyBinIndex] = max(max(results, [], 2));                      % Look for correlation peaks 2D 
    [peakSize,codePhase] = max(max(results));                              % Find code phase of the same correlation peak
    % When searching the secondary peak, it is necessary to exclude the 
    % code piece data around the main peak    
    samplesPerCodeChip = round(fs / CodeRate);                             % how many sample bits in 1 code piece width 
    ExRangeL = codePhase - samplesPerCodeChip;                             % mininum exclude range
    ExRangeR = codePhase + samplesPerCodeChip;                             % maxinum exclude range
    if ExRangeL < 0                                                        % get phase range
        InRange = ExRangeR : BitsPerCode + ExRangeL;                       % if ExRangeL beyonds the index left boundary
    elseif ExRangeR > BitsPerCode                                          % if ExRangeR beyonds the index right boundary
        InRange = ExRangeR - BitsPerCode : ExRangeL;
    else                                                                   % default
        InRange = [1:ExRangeL,ExRangeR : BitsPerCode];
    end
    secondPeakSize = max(results(frequencyBinIndex, InRange));             % Find the second highest correlation peak in the same freq. bin
    acqResults.peakMetric(PRN) = peakSize/secondPeakSize;                  % Store result
    NoiseMean = mean(results(frequencyBinIndex, InRange));                 % calculate noisebase magnitude
    acqResults.PeakNoiseR(PRN) = 10*log10(peakSize/NoiseMean);             % Add a Peak to Noise ratio
    Phase = (1:BitsPerCode)./fs;                                           % mesh the 2D correlation peak
    [TempF,Phase] = meshgrid(1:TempF,Phase);                               % mesh grid
    mesh(TempF,Phase,results');                                            % draw a new axes
    % Fine resolution frequency search
    if (peakSize/secondPeakSize) > settings.acqThreshold                   % If the result is above threshold, then there is a signal ...
        acqResults.acq_Sat_No=acqResults.acq_Sat_No+1;                     % add a new satellite 
        fprintf('%02d ', PRN);                                             % Indicate PRN number of the detected signal 
        caCode = generateCAcode(PRN);                                      % Generate 10msec long C/A codes sequence for given PRN　 
        codeValueIndex = floor((ts * (1:10*BitsPerCode)) / (1/CodeRate));  % The index of the received data aligned at the start position
        longCaCode = caCode((rem(codeValueIndex, 1023) + 1));              % The index ranges from 1 to 1023
        %--- Remove C/A code modulation from the original signal ----------
        % (Using detected C/A code phase)
        xCarrier = signal0DC(codePhase:...
                            (codePhase+10*BitsPerCode-1)).*longCaCode;     % Read 10 ms data
        fftNumPts = 8*(2^(nextpow2(length(xCarrier))));                    % Find the next highest power of two and increase by 8x
        %--- Compute the magnitude of the FFT, find maximum and the
        %associated carrier frequency 
        fftxc = abs(fftshift(fft(xCarrier, fftNumPts))); 
        [~, fftMaxIndex] = max(fftxc);
        fftFreqBins = linspace(-fs/2,fs/2,fftNumPts);
        %--- Save properties of the detected satellite signal -------------
        acqResults.carrFreq(PRN)  = fftFreqBins(fftMaxIndex);
        acqResults.codePhase(PRN) = codePhase;
        acqResults.debugFig = fig;
%         fftxc = abs(fft(xCarrier, fftNumPts)); 
%         
%         uniqFftPts = ceil((fftNumPts + 1) / 2);
%         [~, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5));
%         
%         fftFreqBins = (0 : uniqFftPts-1) * settings.samplingFreq/fftNumPts ;
%         
%         --- Save properties of the detected satellite signal -------------
%         acqResults.carrFreq(PRN)  = fftFreqBins(fftMaxIndex);
%         acqResults.codePhase(PRN) = codePhase;
    else
        fprintf('. ');                                                     % No signal with current PRN 
    end
end 
fprintf(')\n');
fprintf('%d satellites captured',acqResults.acq_Sat_No);
